%% 入射角
%% 计算太阳赤纬角 delta
D=[-59 -28 0 31 61 92 122 153 184 214 245 275 ];
R_total=[];
for i=2:1:size(D,1)
    R=Main(D(i));
    R_total=[R_total;R];
end
fprintf('eta_trunc:%f\n',R);
function R = Main(Day)
D = Day;
STT=[9 10.5 12 13.5 15]; %当地时间,假设为 15 点
result=[];%日结果详细矩阵 5*4 5 个时间点 4 个参数 R 才是最终值
eta_ref=0.92;%镜面反射率
H=80;%吸收塔高度
O_R=[0 0 H];%集热器中心坐标
%定日镜选取根据与中心的距离进行排序
%load 1.mat %经过排序后的定日镜坐标 data---1745*3 x,y,r
data = readtable("附件.xlsx");
data = table2array(data);
h=4;%定日镜安装高度
% 计算太阳赤纬角 delta
sin_delta = sin(2*pi*D/365) * sin(2*pi/360 * 23.45);
delta = asin(sin_delta); % 得到赤纬角的弧度值
phi = deg2rad(39.4); % 北纬 39.4 度，转换为弧度
for XH_ST=1:size(STT,2)
    ST=STT(XH_ST);
    omega =pi/12*(ST-12); % 声明 omega 的值
    % 计算 sin(alpha_s)
    sin_alpha_s = cos(delta) * cos(phi) * cos(omega) + sin(delta) * sin(phi);
    G_0=1.366; %太阳常数
    Hh=3;%3km
    ha=0.4237 - 0.00821*(6 - Hh)^2;%大气压
    hb=0.5055 + 0.00595*(6.5 - Hh)^2;%大气压
    hc=0.2711 + 0.01858*(2.5 - Hh)^2;%大气压
    DNI = G_0*(ha + hb*exp(-hc/sin_alpha_s)); %直接法辐照度
    alpha_s = asin(sin_alpha_s); % 得到 alpha_s 的弧度值
    % 根据公式计算 cos(gamma_s)
    cos_gamma_s = (sin(delta) - sin_alpha_s * sin(phi)) / (cos(alpha_s) * cos(phi));
    gamma_s = acos(cos_gamma_s); % 得到 gamma_s 的弧度值
    x=cos(alpha_s)*sin(gamma_s);
    y=cos(alpha_s)*cos(gamma_s);
    z=sin(alpha_s);
    S_i=[x y z];%太阳入射方向单位向量的相反
    Store=zeros(1745,4);%光学效率 阴影遮挡效率 余弦效率 截断效率
    E_field=0;%定日镜场的输出热功率
    for num=1:1:1745%num=1745;%定日镜编号
        O_A=[data(num,1) data(num,2) h];%定日镜中心坐标
        %计算集热器截断效率 eta_trunc
        d_HR=sqrt((O_A(1)-O_R(1))^2+(O_A(2)-O_R(2))^2+(O_A(3)-O_R(3))^2);%定日镜镜面中心与集热器中心的距离
        eta_at = 0.99321 - 0.0001176*d_HR + 1.97e-8 * (d_HR^2); % 计算吸收塔的吸收率
        S_AR=(O_R-O_A)/norm(O_R-O_A);%定日镜中心到集热器中心的单位向量 %% 反射角
        %% 计算法向量
        S_n=(S_i+S_AR)/norm(S_i+S_AR);%法向量单位向量
        % 计算方位角
        A_H = atan2(S_n(1), S_n(2));
        % 计算俯仰角
        E_H = acos(S_n(3));
        %变换矩阵
        %T=[-sin(E_H),-sin(A_H)*cos(E_H),cos(A_H)*cos(E_H);cos(E_H),-sin(A_H)*sin(E_H),cos(A_H)*sin(E_H);0,cos(A_H),sin(A_H)];
        M1 = [1, 0, 0; 0, cos(E_H), sin(E_H); 0, -sin(E_H), cos(E_H)];
        M2 = [cos(pi-A_H), sin(pi-A_H), 0; -sin(pi-A_H), cos(pi-A_H), 0; 0, 0, 1];
        % 计算矩阵乘积
        T = M1 * M2; %镜坐标系变为地坐标系的转换矩阵
        S_An=[0 0 1];%镜面 A 坐标系上的的法向量
        SSS=S_An*T;
        %计算余弦效率 eta_cos
        eta_cos = dot(S_i, S_n) / (norm(S_i) * norm(S_n)); % 计算入射向量相反和法向量之间的夹角的余弦值
        %计算阴影遮挡效率 eta_sb        
        W_a = 6;
        W_b = 6;
        dx = 0.1;
        dy = 0.1;
        R = 100;
        X = -1/2 * W_a:dx:1/2 * W_a;
        Y = -1/2 * W_b:dy:1/2 * W_b;
        T_inv = inv(T);
        % 在指定的半径内寻找定日镜
        distances = sqrt((data(1:num-1, 1) - O_A(1)).^2 + (data(1:num-1, 2) - O_A(2)).^2);
        valid_indices = find(distances <= R);
        % Main loops
        valid_point_count = 0;%有效点
        valid_points = 0;%每个点的有效光线
        obstructed_mirrors = [];
        for x = X
            for y = Y
                dot_A = [x, y, 0];
                dot_AG = dot_A * T + O_A;
                is_valid = true;
                for idx = valid_indices'
                    O_B = [data(idx, 1), data(idx, 2), h];
                    dot_B = (dot_AG - O_B)*T_inv;
                    S_iB = S_i*T_inv;
                    H2 = [(S_iB(3)*dot_B(1) - S_iB(1)*dot_B(3))/S_iB(3), (S_iB(3)*dot_B(2) -S_iB(2)*dot_B(3))/S_iB(3), 0];
                    if ~(-1/2*W_a <= H2(1) && H2(1) <= 1/2*W_a && -1/2*W_b <= H2(2) && H2(2)<= 1/2*W_b)%不在该平面内为真
                        S_ARB = S_AR*T_inv;
                        H1 = [(S_ARB(3)*dot_B(1) - S_ARB(1)*dot_B(3))/S_ARB(3),(S_ARB(3)*dot_B(2) - S_ARB(2)*dot_B(3))/S_ARB(3), 0];
                        if -1/2*W_a <= H1(1) && H1(1) <= 1/2*W_a && -1/2*W_b <= H1(2) && H1(2)<= 1/2*W_b %在该平面内为真
                            obstructed_mirrors = [obstructed_mirrors; O_B, distances(idx)];
                            is_valid = false;
                            break;
                        end
                    else%在该平面内为真
                        obstructed_mirrors = [obstructed_mirrors; O_B, distances(idx)];
                        is_valid = false;
                        break;
                    end
                end
                if is_valid
                    valid_point_count = valid_point_count + 1;
                    %这个点是有效的计算集热器接收到的光
                    % 计算入射向量和法向量之间的夹角
                    cosTheta_i = dot(S_i, S_n) / (norm(S_i) * norm(S_n));
                    S_nt = S_n / cosTheta_i;
                    d_ki = S_nt - S_AR;
                    d_ki = d_ki / norm(d_ki);
                    theta_range = linspace(-4.65e-3,4.65e-3, 5);
                    theta2_range = linspace(0, 2*pi, 12);
                    % 循环 theta 值
                    for theta = theta_range
                        d_ki2 = tan(theta) * d_ki;
                        % 循环不同的 theta2 值
                        for theta2 = theta2_range                          
                            RR = rotationMatrix(S_AR, -theta2); %顺时针
                            d_tki = RR * d_ki2';
                            d_tki=d_tki';
                            r = d_tki + S_AR;
                            % 解出 t 的范围 判断 r 是否与集热器相交 地系坐标
                            a = r(1)^2 + r(2)^2;
                            b = 2 * (dot_AG(1) * r(1) + dot_AG(2) * r(2));
                            c = dot_AG(1)^2 + dot_AG(2)^2 - 49/4;
                            discriminant = b^2 - 4*a*c;
                            if discriminant >= 0
                                t1 = (-b + sqrt(discriminant)) / (2*a);
                                t2 = (-b - sqrt(discriminant)) / (2*a);
                                % 使用 t 的范围得到 z 的范围 地系坐标                             
                               z1 = dot_AG(3) + t1 * r(3);
                               z2 = dot_AG(3) + t2 * r(3);
                               % 检查 z 的范围是否与(H-4,H+4]有交集
                               z_min = min(z1, z2);
                               z_max = max(z1, z2);
                               if (z_min <= (H+4) && z_max > (H-4))%有交集
                                   valid_points = valid_points + 1;
                               end
                            end
                        end
                    end
                end
            end
        end
        eta_sb = valid_point_count / (length(X) * length(Y));
        Dot_Sum_light=length(theta_range)*length(theta2_range);%每个区域总光线数量
        eta_trunc=valid_points/(valid_point_count*Dot_Sum_light);
        if eta_trunc<1
            disp(num);
            disp(eta_trunc);
        end
        eta=eta_cos*eta_sb*eta_at*eta_trunc*eta_ref;%计算光学效率 eta
        Store(num,1)=eta;
        Store(num,2)=eta_sb;
        Store(num,3)=eta_cos;
        Store(num,4)=eta_trunc;
    end
    %计算定日镜场的输出热功率 E_field
    Ai=36;
    E_field=E_field+Ai*Store(num,1);%E_field=E_field+Ai*eta;
E_field=E_field*DNI;
result=[result;mean(Store) E_field/(num*W_a*W_b)];
end
R=mean(result);
end
function [P_top_1_intersect,P_top_2_intersect,P_bot_1_intersect,P_bot_2_intersect] =shadow_range(i, h) %i 是-S_i
i = i / norm(i);
% 圆柱参数
H = 80; % 吸收塔高度
R = 7/2; % 圆柱半径
O_top= [0 0 (H + R)]; % 圆柱上底面中心坐标
% 计算与入射光 i 垂直的直径方向
diameter_direction = cross(i, [0, 0, 1]);
diameter_direction = diameter_direction / norm(diameter_direction); % 单位化
% 上底面和下底面上与入射光垂直的直径的两个端点
P_top_1 = O_top + R * diameter_direction;
P_top_2 = O_top - R * diameter_direction;
P_bot_1 = P_top_1 - [0, 0, 2*R];
P_bot_2 = P_top_2 - [0, 0, 2*R];
% 计算在 z=h 时的交点
t_top_1 = (h - P_top_1(3)) / i(3);
P_top_1_intersect = P_top_1 + t_top_1 * i;
t_top_2 = (h - P_top_2(3)) / i(3);
P_top_2_intersect = P_top_2 + t_top_2 * i;
t_bot_1 = (h - P_bot_1(3)) / i(3);
P_bot_1_intersect = P_bot_1 + t_bot_1 * i;
t_bot_2 = (h - P_bot_2(3)) / i(3);
P_bot_2_intersect = P_bot_2 + t_bot_2 * i;
end
function isInside = pointInRectangle(x, y, P_top_1_intersect, P_top_2_intersect,P_bot_1_intersect, P_bot_2_intersect)
% 定义矩形的顶点
polygon = [P_top_1_intersect; P_top_2_intersect; P_bot_2_intersect; P_bot_1_intersect;P_top_1_intersect];
% 计算交点数
crossings = 0;
for i = 1:length(polygon)-1
    if ((polygon(i, 2) > y) ~=(polygon(i+1, 2) > y))&&(x < (polygon(i+1, 1) - polygon(i, 1)) * (y - polygon(i, 2)) /(polygon(i+1, 2) - polygon(i, 2)) + polygon(i, 1))
        crossings = crossings + 1;
    end
end
% 根据交点数判断点是否在矩形内
if mod(crossings, 2) == 1
    isInside = true;
else
    isInside = false;
end
function RR = rotationMatrix(axis, angle)
c = cos(angle);
s = sin(angle);
t = 1 - c;
x = axis(1);
y = axis(2);
z = axis(3);
RR = [t*x*x + c, t*x*y - s*z, t*x*z + s*y;t*x*y + s*z, t*y*y + c, t*y*z - s*x;t*x*z - s*y, t*y*z + s*x, t*z*z + c];
end
end